
/*------------------------------------------------------------------------
 *  Copyright (C) 2010  Luis M. de la Cruz
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _SUBDOMAIN_HPP_
#define _SUBDOMAIN_HPP_


#include "DDM/CartComm.hpp"
#include "Meshes/StructuredMesh.hpp"

namespace Tuna {
  
template<typename Tprec, int Dim>
class SubDomain {
 
public: 
  typedef Tprec prec_t;
  typedef typename TunaArray<prec_t, Dim>::huge ScalarField;
  typedef typename TunaArray<prec_t,1>::huge array1D;

  SubDomain(CartComm<Dim>& cart) {
    neighbor    = cart.getNeighbors(); 
    I           = cart.get_I();
    J           = cart.get_J();
    K           = cart.get_K();
    rank        = cart.get_rank();
    nsubdoms[0] = cart.getNumProc_I();
    nsubdoms[1] = cart.getNumProc_J();
    nsubdoms[2] = cart.getNumProc_K();

    cart_p = &cart;

    for(int i = 0; i < 6; i++) {
      overlap[i] = 0;
      cells_overlap[i] = 1;
      send_indx[i] = 0;
    }

  }

  SubDomain(CartComm<Dim>& cart, StructuredMesh<Uniform<double, Dim> >& mesh,
	    int nc_ovlp_l, int nc_ovlp_r,
	    int nc_ovlp_d = 0, int nc_ovlp_u = 0,
	    int nc_ovlp_b = 0, int nc_ovlp_f = 0) {
    neighbor    = cart.getNeighbors(); 
    I           = cart.get_I();
    J           = cart.get_J();
    K           = cart.get_K();
    rank        = cart.get_rank();
    nsubdoms[0] = cart.getNumProc_I();
    nsubdoms[1] = cart.getNumProc_J();
    nsubdoms[2] = cart.getNumProc_K();

    cart_p = &cart;

    for(int i = 0; i < 6; i++) {
      overlap[i] = 0;
      cells_overlap[i] = 1;
      send_indx[i] = 0;
    }
 
// Global Mesh
//  0     1     2     3     4     5     6     7    8      9     10    11   12
//  |                                   |                                   |
//  +--*--+--*--+--*--+--*--+--*--+--*--+--*--+--*--+--*--+--*--+--*--+--*--+
//  |                                   |                                   |
//  0  1     2     3     4     5     6     7     8     9    10    11    12 13
//   
// Local mesh for subdomains of equal length    
//  0     1     2     3     4     5     6     7
//  |                                   |     
//  +--*--+--*--+--*--+--*--+--*--+--*--+--*--+
//  |                                   |
//  0  1     2     3     4     5     6     7
//
//                                0     1     2     3     4     5     6     7
//                                      |                                   |
//                                +--*--+--*--+--*--+--*--+--*--+--*--+--*--+
//                                      |                                   |
//                                   0     1     2     3     4     5     6  7  
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    double dz = mesh.getDelta(Z);

    double length_x = mesh.getLength(X) / nsubdoms[0] // Non-overlapping length
      + createOverlap(LEFT,  nc_ovlp_l, dx)    // Overlap on the LEFT
      + createOverlap(RIGHT, nc_ovlp_r, dx);   // Overlap on the RIGHT
    
    double length_y = mesh.getLength(Y) / nsubdoms[1] // Non-overlapping length
      + createOverlap(DOWN, nc_ovlp_d, dy)    // Overlap on the DOWN
      + createOverlap(UP,   nc_ovlp_u, dy);   // Overlap on the UP

    double length_z = mesh.getLength(Z) / nsubdoms[2] // Non-overlapping length
      + createOverlap(BACK, nc_ovlp_b, dz)    // Overlap on the BACK
      + createOverlap(FRONT, nc_ovlp_f, dz);   // Overlap on the FRONT

    // Restriction: (num_nodes_? - 1) % num_subdom = 0 
    int num_nodes_x = (mesh.getExtentNodes(X) - 1) / nsubdoms[0] + 1 // Non-overlapping part
      + nc_ovlp_l  // Number of overlapping nodes on the left  
      + nc_ovlp_r; // Number of overlapping nodes on the right  

    int num_nodes_y = (mesh.getExtentNodes(Y) - 1) / nsubdoms[1] + 1 // Non-overlapping part
      + nc_ovlp_d  // Number of overlapping nodes on the down  
      + nc_ovlp_u; // Number of overlapping nodes on the up 

    int num_nodes_z = (mesh.getExtentNodes(Z) - 1) / nsubdoms[2] + 1 // Non-overlapping part
      + nc_ovlp_b  // Number of overlapping nodes on the back
      + nc_ovlp_f; // Number of overlapping nodes on the front 

    // Resizing the mesh to the local values for the subdomain
    if (Dim == 1)
      mesh.resize(length_x, num_nodes_x);
    if (Dim == 2)
      mesh.resize(length_x,  num_nodes_x, length_y, num_nodes_y);
    if (Dim == 3)
      mesh.resize(length_x,  num_nodes_x, length_y, num_nodes_y, length_z, num_nodes_z);

    // Define the info to be exchanged between subdomains
    setupInfoExchange(mesh);
  }

  ~SubDomain() { };

  inline prec_t createOverlap(cartDir_t n, int& ov, double delta) { 

    cells_overlap[n] = ov;
    overlap[n] = delta * ov;

    // Check for the boundary subdomains...
    if( 
       ( (I == 0)               && (n == LEFT)  ) ||
       ( (I == nsubdoms[0] - 1) && (n == RIGHT) ) ||
       ( (J == 0)               && (n == DOWN)  ) ||
       ( (J == nsubdoms[1] - 1) && (n == UP)    ) || 
       ( (K == 0)               && (n == BACK)  ) ||
       ( (K == nsubdoms[2] - 1) && (n == FRONT) )   
      )     
      {    
	overlap[n]       = 0;
	cells_overlap[n] = 0;
	ov               = 0;
      }
    
    return overlap[n];
  }

  /*!
   *
   * \verbatim
   
   |                     Original Domain                       |
   +----+----+----+----+----+----+----+----+----+----+----+----+
   |                                                           |

                        Domain Decomposition
   |         L                                        R         |
   +----+----+----+----+....+         +....+----+----+----+----+
                          ^             ^
                          |             |
                     c_ovlp_r_L     c_ovlp_l_R


                       |         C         |     
                  +....+----+----+----+----+....+
                     ^                       ^
                     |                       |
            cells_overlap[LEFT]      cells_overlap[RIGHT]

     \endverbatim
 * 
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Sat May 14 11:28:46 CDT 2011 ]
 ***************************************************************************
 */ 

  template<typename Tmesh>
  inline bool setupInfoExchange(Tmesh& mesh)
  {
    int Nx = mesh.getExtentNodes(X); //< Number of nodes on x-axis
    int Ny = mesh.getExtentNodes(Y); //< Number of nodes on y-axis
    int Nz = mesh.getExtentNodes(Z); //< Number of nodes on z-axis

    int c_ovlp_r_L = 0; //< cells on the right overlap of LEFT neighbor
    int c_ovlp_l_R = 0; //< cells on the left overlap of RIGHT neighbor

    // Exchange overlapping in x-axis direction 
    /**
    cart_p->comm.Sendrecv(&c_ovlp_r, 1, MPI::INT,  // Send to right
			  neighbor[RIGHT], 0,
			  &c_ovlp_r_l, 1, MPI::INT, // Receive from left
			  neighbor[LEFT], 0);
    **/
    MPI_Sendrecv(&cells_overlap[RIGHT], 1, MPI::INT, neighbor[RIGHT], 0, // Send to right
		 &c_ovlp_r_L, 1, MPI::INT, neighbor[LEFT], 0, // Receive from left
		 cart_p->comm, &status);
    /**
    cart_p->comm.Sendrecv(&c_ovlp_l, 1, MPI::INT, // Send to left 
			  neighbor[LEFT], 0,
			  &c_ovlp_l_r, 1, MPI::INT, // Receive from right
			  neighbor[RIGHT], 0);
    **/

    MPI_Sendrecv(&cells_overlap[LEFT], 1, MPI::INT, neighbor[LEFT], 0, // Send to left 
		 &c_ovlp_l_R, 1, MPI::INT, neighbor[RIGHT], 0, // Receive from right
		 cart_p->comm, &status);

    send_indx[LEFT] = c_ovlp_r_L + cells_overlap[LEFT];
    send_indx[RIGHT] = Nx - 1 - ( c_ovlp_l_R + cells_overlap[RIGHT] );

    /* *
    std::cout << "\n I = " << I << " s_l = " << send_indx[LEFT]
	      << " s_r = " << send_indx[RIGHT] << "\n\n";
    /* */    

    int c_ovlp_u_D = 0; //< cells on the up overlap of DOWN neighbor
    int c_ovlp_d_U = 0; //< cells on the down overlap of UP neighbor
    
    // Exchange overlapping in y-axis direction
    if( Dim >= 2) {
      MPI_Sendrecv(&cells_overlap[UP], 1, MPI::INT, neighbor[UP], 0, // Send to up
		   &c_ovlp_u_D, 1, MPI::INT, neighbor[DOWN], 0, // Receive from down
		   cart_p->comm, &status);

      MPI_Sendrecv(&cells_overlap[DOWN], 1, MPI::INT, neighbor[DOWN], 0, // Send to down 
		   &c_ovlp_d_U, 1, MPI::INT, neighbor[UP], 0, // Receive from up
		   cart_p->comm, &status);
      
      send_indx[DOWN] = c_ovlp_u_D + cells_overlap[DOWN];
      send_indx[UP] = Ny - 1 - ( c_ovlp_d_U + cells_overlap[UP] );  

    /* *
    std::cout << "\n J = " << J << " s_d = " << send_indx[DOWN]
	      << " s_u = " << send_indx[UP] << "\n\n";
    /* */
    }

    int c_ovlp_f_B = 0; //< cells on the front overlap of BACK neighbor
    int c_ovlp_b_F = 0; //< cells on the back overlap of FRONT neighbor
    
    // Exchange overlapping in y-axis direction
    if( Dim >= 3) {
      MPI_Sendrecv(&cells_overlap[FRONT], 1, MPI::INT, neighbor[FRONT], 0, // Send to FRONT
		   &c_ovlp_f_B, 1, MPI::INT, neighbor[BACK], 0, // Receive from BACK
		   cart_p->comm, &status);

      MPI_Sendrecv(&cells_overlap[BACK], 1, MPI::INT, neighbor[BACK], 0, // Send to BACK 
		   &c_ovlp_b_F, 1, MPI::INT, neighbor[FRONT], 0, // Receive from FRONT
		   cart_p->comm, &status);
      
      send_indx[BACK] = c_ovlp_f_B + cells_overlap[BACK];
      send_indx[FRONT] = Nz - 1 - ( c_ovlp_b_F + cells_overlap[FRONT] );  

    /* *
    std::cout << "\n K = " << K << " s_f = " << send_indx[FRONT]
	      << " s_b = " << send_indx[BACK] << "\n\n";
    /* */
    }


    return 0;
  }

  template<typename T_eq>
  inline void setRealDirichlet(T_eq& eq, BC_t wall, prec_t w_v = 0.0);

  template<typename T_eq>
  inline void setRealNeumann(T_eq& eq, BC_t wall, prec_t w_v = 0.0);

  template<typename T_eq>
  inline void setGhostDirichlet(T_eq& eq);

  template<typename T_eq>
  inline void setGhostNeumann(T_eq& eq, prec_t w_v = 0.0);

  template<typename T_eq>
  inline void setGhost(T_eq& eq, prec_t w_v = 0.0);

  inline bool infoExchange(ScalarField& phi, double alpha);
  inline bool infoExchange1(ScalarField& phi);


  inline prec_t getOverlap(cartDir_t n) { return overlap[n]; }

  void print () {
    std::cout << "\n +----- TUNA info -----+"
	      << "\n | SubDomain < " << Dim << "  > "
	      << "\n +-----+";
    
    if (Dim == 1) {
      std::cout << "\n | Proc " << rank << " (" << I << ") "
		<< " | RIGHT = " << neighbor[RIGHT]
		<< " , " << cells_overlap[RIGHT]
		<< " , " << overlap[RIGHT]
		<< " | LEFT = " << neighbor[LEFT]
		<< " , " << cells_overlap[LEFT]
		<< " , " << overlap[LEFT];
      }

      if (Dim == 2) {
	std::cout << "\n | Proc " << rank << " (" << I << "," << J << ") "
		  << " | RIGHT = " << neighbor[RIGHT]
		  << " , " << cells_overlap[RIGHT]
		  << " , " << overlap[RIGHT]
		  << " | LEFT = " << neighbor[LEFT]
		  << " , " << cells_overlap[LEFT]
		  << " , " << overlap[LEFT]
		  << " | UP = " << neighbor[UP]
		  << " , " << cells_overlap[UP]
		  << " , " << overlap[UP]
		  << " | DOWN = " << neighbor[DOWN]
		  << " , " << cells_overlap[DOWN]
		  << " , " << overlap[DOWN];       
      }

      if (Dim == 3) {
	std::cout << "\n | Proc " << rank << " (" << I << "," << J << "," << K << ") "
		  << " | RIGHT = " << neighbor[RIGHT]
		  << " , " << cells_overlap[RIGHT]
		  << " , " << overlap[RIGHT]
		  << " | LEFT = " << neighbor[LEFT]
		  << " , " << cells_overlap[LEFT]
		  << " , " << overlap[LEFT]
		  << " | UP = " << neighbor[UP]
		  << " , " << cells_overlap[UP]
		  << " , " << overlap[UP]
		  << " | DOWN = " << neighbor[DOWN]
		  << " , " << cells_overlap[DOWN]
		  << " , " << overlap[DOWN]
		  << " | FRONT = " << neighbor[FRONT]
		  << " , " << cells_overlap[FRONT]
		  << " , " << overlap[FRONT]
		  << " | BACK = " << neighbor[BACK]
		  << " , " << cells_overlap[BACK]
		  << " , " << overlap[BACK];
      }

      std::cout << "\n +---------------------+\n "; 
    }

private:
  prec_t overlap[6];    // Size of the overlap on each side
  int cells_overlap[6]; // number of cells that overlap on each side
  int send_indx[6];     // Indexes of cells to send on each side
  const int* neighbor;  // Neighbors identification
  int I, J, K, rank, nsubdoms[3];
  CartComm<Dim> *cart_p;
  MPI_Status status;
};
  


template<typename Tprec, int Dim>
inline bool SubDomain<Tprec, Dim>::infoExchange1(ScalarField& phi)
{

  if (Dim == 1) {
    if ( nsubdoms[0] > 1) {  
	double phi_send_right = ( phi(send_indx[RIGHT]) + 
				  phi(send_indx[RIGHT] + 1) ) * 0.5;
	double phi_send_left =  ( phi(send_indx[LEFT]) + 
				  phi(send_indx[LEFT] + 1) ) * 0.5;

	double *ri_send, *le_send, *ri_recv, *le_recv;
	ri_send = &phi_send_right;
	le_send = &phi_send_left;
	ri_recv = &phi( phi.ubound(firstDim) );
	le_recv = &phi( phi.lbound(firstDim) );

	cart_p->comm.Sendrecv(ri_send, 1, MPI::DOUBLE, // Send info to the right
			      neighbor[RIGHT], 0,
			      le_recv, 1, MPI::DOUBLE, // Receive info from
			      neighbor[LEFT], 0);   // the left
	
	cart_p->comm.Sendrecv(le_send, 1, MPI::DOUBLE, // Send info to the right
			      neighbor[LEFT], 0,
			      ri_recv, 1, MPI::DOUBLE, // Receive info from 
			      neighbor[RIGHT], 0); // the left
	

	if ( I != (nsubdoms[0] - 1) )
	  phi( phi.ubound(firstDim) ) = *ri_recv;
	if (I != 0) 
	  phi( phi.lbound(firstDim) ) = *le_recv;
    }      
  }

  if (Dim == 2) {
    
    int Nx = phi.extent(firstDim);
    int Ny = phi.extent(secondDim);

    MPI::Datatype Vector_x = MPI::DOUBLE.Create_contiguous(Nx);
    MPI::Datatype Vector_y = MPI::DOUBLE.Create_contiguous(Ny);    
    Vector_x.Commit();
    Vector_y.Commit();

    if ( nsubdoms[0] > 1) {  
      
      double ri_send[Ny], le_send[Ny], ri_recv[Ny], le_recv[Ny];
      
      if (I != 0)     
	for(int j = 0; j < Ny; j++) 
	  le_send[j] = ( phi(send_indx[LEFT]    , j) +
			 phi(send_indx[LEFT] + 1, j) ) * 0.5;
      
      if (I != (nsubdoms[0] - 1) ) 
	for(int j = 0; j < Ny; j++) 
	  ri_send[j] = ( phi(send_indx[RIGHT]    , j) +
			 phi(send_indx[RIGHT] + 1, j) ) * 0.5;
      
      cart_p->comm.Sendrecv(le_send, 1, Vector_y, neighbor[LEFT], 0,
			    ri_recv, 1, Vector_y, neighbor[RIGHT], 0);
      
      cart_p->comm.Sendrecv(ri_send, 1, Vector_y, neighbor[RIGHT], 0,
			    le_recv, 1, Vector_y, neighbor[LEFT], 0);
      
      if (I != 0) 
	for(int j = 0; j < Ny; j++) 
	  phi(0    ,j) = le_recv[j];
      
      if (I != (nsubdoms[0] - 1) ) 
	for(int j = 0; j < Ny; j++) 
	  phi(Nx-1 ,j) = ri_recv[j];
    }
    
    if ( nsubdoms[1] > 1) {
      
      double up_send[Nx], do_send[Nx], up_recv[Nx], do_recv[Nx];
      
      if (J != 0)     
	for(int i = 0; i < Nx; i++) 
	  do_send[i] = ( phi(i, send_indx[DOWN]) +
			 phi(i, send_indx[DOWN] + 1) ) * 0.5;
      
      if (J != (nsubdoms[1] - 1) ) 
	for(int i = 0; i < Nx; i++) 
	  up_send[i] = ( phi(i, send_indx[UP]) + 
			 phi(i, send_indx[UP] + 1) ) * 0.5;
      
      
      cart_p->comm.Sendrecv(do_send, 1, Vector_x, neighbor[DOWN], 1,
			    up_recv, 1, Vector_x, neighbor[UP], 1);
      
      cart_p->comm.Sendrecv(up_send, 1, Vector_x, neighbor[UP], 1,
			    do_recv, 1, Vector_x, neighbor[DOWN], 1);
            
      if (J != 0) 
	for(int i = 0; i < Nx; i++) 
	  phi(i, 0) = do_recv[i];
      
      if (J != (nsubdoms[1] - 1) ) 
	for(int i = 0; i < Nx; i++) 
	  phi(i, Ny-1) = up_recv[i];
      
    }

    Vector_x.Free();
    Vector_y.Free();
  }

  if (Dim == 3) {
  }
    return 0;
}

template<typename Tprec, int Dim>
inline bool SubDomain<Tprec, Dim>::infoExchange(ScalarField& phi, double alpha = 0.5)
{

  //
  // alpha = 0 you get the value on nodes.
  // alpha = 0.5 is for FVM on uniform meshes.
  // alpha = \Delta x_i / (\Delta x_i + \Delta x_i+1) for non-uniform meshes
  //

  if (Dim == 1) {
    if ( nsubdoms[0] > 1) { 

      double alphaR = alpha, alphaL = alpha;
      double phi_send_right = phi(send_indx[RIGHT]    ) * (1 - alphaR) + 
	                      phi(send_indx[RIGHT] + 1) * alphaR;
      double phi_send_left =  phi(send_indx[LEFT]    ) * (1 - alphaL)+ 
			      phi(send_indx[LEFT] + 1) * alphaL;
 
	double *ri_send, *le_send, *ri_recv, *le_recv;
	ri_send = &phi_send_right; // &phi( send_indx[RIGHT] );
	le_send = &phi_send_left;  // &phi( send_indx[LEFT] );
	ri_recv = &phi( phi.ubound(firstDim) );
	le_recv = &phi( phi.lbound(firstDim) );

	/**
	cart_p->comm.Sendrecv(ri_send, 1, MPI::DOUBLE, // Send info to the right
			      neighbor[RIGHT], 0,
			      le_recv, 1, MPI::DOUBLE, // Receive info from
			      neighbor[LEFT], 0);   // the left
	**/
	MPI_Sendrecv(ri_send, 1, MPI::DOUBLE, neighbor[RIGHT], 0, // Send info to the right
		     le_recv, 1, MPI::DOUBLE, neighbor[LEFT], 0, // Recv info from the left
		     cart_p->comm, &status);

	MPI_Sendrecv(le_send, 1, MPI::DOUBLE, neighbor[LEFT], 0, // Send info to the right
		     ri_recv, 1, MPI::DOUBLE, neighbor[RIGHT],0, // Recv info from the left
		     cart_p->comm, &status);

	if ( I != (nsubdoms[0] - 1) )
	  phi( phi.ubound(firstDim) ) = *ri_recv;
	if (I != 0) 
	  phi( phi.lbound(firstDim) ) = *le_recv;
	
    }      
  }


  if (Dim == 2) {

    int Nx = phi.extent(firstDim);
    int Ny = phi.extent(secondDim);

    MPI::Datatype Vector_x = MPI::DOUBLE.Create_contiguous(Nx);
    MPI::Datatype Vector_y = MPI::DOUBLE.Create_contiguous(Ny);    
    Vector_x.Commit();
    Vector_y.Commit();
    
    if ( nsubdoms[0] > 1) {  
      
      double alphaL = alpha, alphaR = alpha;
      double ri_send[Ny], le_send[Ny], ri_recv[Ny], le_recv[Ny];
      
      if (I != 0)     
	for(int j = 0; j < Ny; j++)
	  le_send[j] = phi(send_indx[LEFT], j) * (1 - alphaL) +
	               phi(send_indx[LEFT] + 1, j) * alphaL;

      if (I != (nsubdoms[0] - 1) ) 
	for(int j = 0; j < Ny; j++) 
	  ri_send[j] = phi(send_indx[RIGHT], j) * (1 - alphaR)+
		       phi(send_indx[RIGHT] + 1, j) * alphaR;
      

      MPI_Sendrecv(le_send, 1, Vector_y, neighbor[LEFT], 0, // send to left
		   ri_recv, 1, Vector_y, neighbor[RIGHT], 0,  // recv from right
		   cart_p->comm, &status);
      
      MPI_Sendrecv(ri_send, 1, Vector_y, neighbor[RIGHT], 0, // send to right
		   le_recv, 1, Vector_y, neighbor[LEFT], 0, // recv from left
      		   cart_p->comm, &status);

      if (I != 0) 
	for(int j = 0; j < Ny; j++) 
	  phi(0    ,j) = le_recv[j];
      
      if (I != (nsubdoms[0] - 1) ) 
	for(int j = 0; j < Ny; j++) 
	  phi(Nx-1 ,j) = ri_recv[j];
    }
    
    if ( nsubdoms[1] > 1) {
      
      double alphaD = alpha, alphaU = alpha;
      double up_send[Nx], do_send[Nx], up_recv[Nx], do_recv[Nx];
      
      if (J != 0)     
	for(int i = 0; i < Nx; i++) 
	  do_send[i] = phi(i, send_indx[DOWN]) * (1 - alphaD) +
      		       phi(i, send_indx[DOWN] + 1) * alphaD;

      if (J != (nsubdoms[1] - 1) ) 
	for(int i = 0; i < Nx; i++) 
	  up_send[i] = phi(i, send_indx[UP]) * (1 - alphaU) +
                       phi(i, send_indx[UP] + 1) * alphaU;
              

      
      MPI_Sendrecv(do_send, 1, Vector_x, neighbor[DOWN], 1, // send to down
		   up_recv, 1, Vector_x, neighbor[UP], 1,   // recv from up
		   cart_p->comm, &status);
      
      MPI_Sendrecv(up_send, 1, Vector_x, neighbor[UP], 1, // send to up
		   do_recv, 1, Vector_x, neighbor[DOWN], 1, // recv from down
      		   cart_p->comm, &status);
      
      if (J != 0) 
	for(int i = 0; i < Nx; i++) 
	  phi(i, 0) = do_recv[i];
      
      if (J != (nsubdoms[1] - 1) ) 
	for(int i = 0; i < Nx; i++) 
	  phi(i, Ny-1) = up_recv[i];
      
    }

    Vector_x.Free();
    Vector_y.Free();
  }

  /****
  if (Dim == 3) {

    int Nx = phi.extent(firstDim);
    int Ny = phi.extent(secondDim);
    int Nz = phi.extent(thirdDim);

    MPI::Datatype Vector_x = MPI::DOUBLE.Create_contiguous(Ny*Nz);
    MPI::Datatype Vector_y = MPI::DOUBLE.Create_contiguous(Nx*Nz);    
    MPI::Datatype Vector_z = MPI::DOUBLE.Create_contiguous(Nz*Ny);    
    Vector_x.Commit();
    Vector_y.Commit();
    Vector_z.Commit();
    
    if ( nsubdoms[0] > 1) {  
      
      double ri_send[Ny], le_send[Ny], ri_recv[Ny], le_recv[Ny];
      
      if (I != 0)     
	for(int j = 0; j < Ny; j++) 
	  le_send[j] = phi(send_indx[LEFT], j);
      
      if (I != (nsubdoms[0] - 1) ) 
	for(int j = 0; j < Ny; j++) 
	  ri_send[j] = phi(send_indx[RIGHT], j);
      
      MPI_Sendrecv(le_send, 1, Vector_y, neighbor[LEFT], 0, // send to left
		   ri_recv, 1, Vector_y, neighbor[RIGHT], 0,  // recv from right
		   cart_p->comm, &status);
      
      MPI_Sendrecv(ri_send, 1, Vector_y, neighbor[RIGHT], 0, // send to right
		   le_recv, 1, Vector_y, neighbor[LEFT], 0, // recv from left
      		   cart_p->comm, &status);

      if (I != 0) 
	for(int j = 0; j < Ny; j++) 
	  phi(0    ,j) = le_recv[j];
      
      if (I != (nsubdoms[0] - 1) ) 
	for(int j = 0; j < Ny; j++) 
	  phi(Nx-1 ,j) = ri_recv[j];
    }
    
    if ( nsubdoms[1] > 1) {
      
      double up_send[Nx], do_send[Nx], up_recv[Nx], do_recv[Nx];
      
      if (J != 0)     
	for(int i = 0; i < Nx; i++) 
	  do_send[i] = phi(i, send_indx[DOWN]);
      
      if (J != (nsubdoms[1] - 1) ) 
	for(int i = 0; i < Nx; i++) 
	  up_send[i] = phi(i, send_indx[UP]);

      
      MPI_Sendrecv(do_send, 1, Vector_x, neighbor[DOWN], 1, // send to down
		   up_recv, 1, Vector_x, neighbor[UP], 1,   // recv from up
		   cart_p->comm, &status);
      
      MPI_Sendrecv(up_send, 1, Vector_x, neighbor[UP], 1, // send to up
		   do_recv, 1, Vector_x, neighbor[DOWN], 1, // recv from down
      		   cart_p->comm, &status);
      
      if (J != 0) 
	for(int i = 0; i < Nx; i++) 
	  phi(i, 0) = do_recv[i];
      
      if (J != (nsubdoms[1] - 1) ) 
	for(int i = 0; i < Nx; i++) 
	  phi(i, Ny-1) = up_recv[i];
      
    }

    Vector_x.Free();
    Vector_y.Free();
  }
  ***/

  return 0;
}


template<typename Tprec, int Dim>
template<typename T_eq>
inline void SubDomain<Tprec, Dim>::setRealDirichlet(T_eq& eq, 
						    BC_t wall, 
						    prec_t value)
{
  if (I == 0 && wall == LEFT_WALL) {
    eq.setDirichlet(wall, value);
  } else if (I == nsubdoms[0] - 1 && wall == RIGHT_WALL) {
    eq.setDirichlet(wall, value);
  } 

  if (Dim >= 2) {
    if (J == 0 && wall == BOTTOM_WALL) {
      eq.setDirichlet(wall, value);
    } else if (J == nsubdoms[1] - 1 && wall == TOP_WALL) {
      eq.setDirichlet(wall, value);
    } 
  } 

  if (Dim == 3) {
    if (K == 0 && wall == BACK_WALL) {
      eq.setDirichlet(wall, value);
    } else if (K == nsubdoms[2] - 1 && wall == FRONT_WALL) {
      eq.setDirichlet(wall, value);
    } 
  }
}
      
template<typename Tprec, int Dim>
template<typename T_eq>
inline void SubDomain<Tprec, Dim>::setRealNeumann(T_eq& eq,
						  BC_t wall, 
						  prec_t value) 
{
  if (I == 0 && wall == LEFT_WALL) {
    eq.setNeumann(wall, value);
  } else if (I == nsubdoms[0] - 1 && wall == RIGHT_WALL) {
    eq.setNeumann(wall, value);
  } 
  
  if (Dim >= 2) {
    if (J == 0 && wall == BOTTOM_WALL) {
      eq.setNeumann(wall, value);
    } else if (J == nsubdoms[1] - 1 && wall == TOP_WALL) {
      eq.setNeumann(wall, value);
    } 
  } 
  
  if (Dim == 3) {
    if (K == 0 && wall == BACK_WALL) {
      eq.setNeumann(wall, value);
    } else if (K == nsubdoms[2] - 1 && wall == FRONT_WALL) {
      eq.setNeumann(wall, value);
    } 
  }

}

template<typename Tprec, int Dim>
template<typename T_eq>
inline void SubDomain<Tprec, Dim>::setGhostDirichlet(T_eq& eq)
{
  if ( nsubdoms[0] > 1) {
    if (I == 0) {
      eq.setDirichlet(RIGHT_WALL);
    } else if (I == nsubdoms[0] - 1) {
      eq.setDirichlet(LEFT_WALL);
    } else {
      eq.setDirichlet(RIGHT_WALL);
      eq.setDirichlet(LEFT_WALL);
    } 
  }

  if (Dim >= 2) {
    if ( nsubdoms[1] > 1) {
      if (J == 0) {
	eq.setDirichlet(TOP_WALL);
      } else if (J == nsubdoms[1] - 1) {
	eq.setDirichlet(BOTTOM_WALL);
      } else {
	eq.setDirichlet(BOTTOM_WALL);
	eq.setDirichlet(TOP_WALL);
      }
    }
  } 
  
  if (Dim == 3) {
    if ( nsubdoms[2] > 1) {
      if (K == 0) {
	eq.setDirichlet(FRONT_WALL);
      } else if (K == nsubdoms[2] - 1) {
	eq.setDirichlet(BACK_WALL);
      } else {
	eq.setDirichlet(BACK_WALL);
	eq.setDirichlet(FRONT_WALL);
      }
    }
  }
}
    
template<typename Tprec, int Dim>
template<typename T_eq>
inline void SubDomain<Tprec, Dim>::setGhostNeumann(T_eq& eq, prec_t value)
{
  if ( nsubdoms[0] > 1) {
    if (I == 0) {
      eq.setNeumann(RIGHT_WALL, value);
    } else if (I == nsubdoms[0] - 1) {
      eq.setNeumann(LEFT_WALL, value);
    } else {
      eq.setNeumann(RIGHT_WALL, value);
      eq.setNeumann(LEFT_WALL, value);
    } 
  }

  if (Dim >= 2) {
    if ( nsubdoms[1] > 1) {
      if (J == 0) {
	eq.setNeumann(TOP_WALL, value);
      } else if (J == nsubdoms[1] - 1) {
	eq.setNeumann(BOTTOM_WALL, value);
      } else {
	eq.setNeumann(BOTTOM_WALL, value);
	eq.setNeumann(TOP_WALL, value);
      }
    }
  } 
  
  if (Dim == 3) {
    if ( nsubdoms[2] > 1) {
      if (K == 0) {
	eq.setNeumann(FRONT_WALL, value);
      } else if (K == nsubdoms[2] - 1) {
	eq.setNeumann(BACK_WALL, value);
      } else {
	eq.setNeumann(BACK_WALL, value);
	eq.setNeumann(FRONT_WALL, value);
      }
    }
  } 
}


template<typename Tprec, int Dim>
template<typename T_eq>
inline void SubDomain<Tprec, Dim>::setGhost(T_eq& eq, prec_t value)
{
  if ( nsubdoms[0] > 1) {
    if (I == 0) {
      eq.set_ei(eq.get_ei() - 1);
    } else if (I == nsubdoms[0] - 1) {
      eq.set_bi(eq.get_bi() + 1);
    } else {
      eq.set_ei(eq.get_ei() - 1);
      eq.set_bi(eq.get_bi() + 1);
    } 
  }

  if (Dim >= 2) {
    if ( nsubdoms[1] > 1) {
      if (J == 0) {
	eq.set_ej(eq.get_ej() - 1);
      } else if (J == nsubdoms[1] - 1) {
	eq.set_bj(eq.get_bj() + 1);
      } else {
	eq.set_ej(eq.get_ej() - 1);
	eq.set_bj(eq.get_bj() + 1);
      }
    }
  } 
  
  if (Dim == 3) {
    if ( nsubdoms[2] > 1) {
      if (K == 0) {
	eq.set_ek(eq.get_ek() - 1);
      } else if (K == nsubdoms[2] - 1) {
	eq.set_bk(eq.get_bk() + 1);
      } else {
	eq.set_ek(eq.get_ek() - 1);
	eq.set_bk(eq.get_bk() + 1);
      }
    }
  } 
}

} // Tuna namespace

#endif //_SUBDOMAIN_HPP_
